arXiv: 1502.0053 lv2 [physics.plasm-ph] 15 Sep 2015 


Mon. Not. R. Astron. Soc. 000, 1-?? (2015) Printed 16 September 2015 (MN IATeK style file v2.2) 


Particle acceleration and radiation friction effects in the 
filamentation instability of pair plasmas 

M. D’Angelo^*, L. Fedeli^’^, A. Sgattoni^, F. Pegoraro^ and A. Macchi^’^ 

^ Gran Sasso Science Institute-INFN, viale Francesco Crispi 7, L’Aquila, 67100, Italy 
Dipartimento di Fisica Enrico Fermi, Universitd di Pisa, Largo Bruno Pontecorvo 3, Pisa 1-56127, Italy 
^Istituto Nazionale di Ottica, Consiglio Nazionale delle Ricerche (CNR/INO), u.o.s. Adriano Gozzini, Pisa, Italy 


16 September 2015 


ABSTRACT 

The evolution of the filamentation instability produced by two counter-streaming, 
ultrarelativistic pair plasmas is studied with particle-in-cell simulations. Radiation 
friction effects are taken into account. Two dimensional simulations are performed for 
both cases of the initial momenta being perpendicular (T-mode) or parallel (P-mode) 
to the simulation plane. In the initial stage the instability is purely transverse for both 
modes and generates small-scale filaments which later merge into larger structures. 
Particle acceleration leads to a strong broadening of the energy spectrum with the 
formation of a peak at twice the initial energy for the T-mode. In the nonlinear stage 
significant differences between T- and P-modes in the evolution of the fields and in 
the spectra of accelerated particles are apparent. The presence of radiative losses does 
not change the dynamics of the instability but strongly affects the structure of the 
particle spectra in the ultra-relativistic regime (particle energy > 100 MeV) and for 
high plasma densities (> 10^^ cm“^). 
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1 INTRODUCTION 

From the 1970s on, the long-standing problem of high- 
energy cosmic ray origin has involved beam-plasma 
instabilities in order to explain some aspects of the accel¬ 
eration mechanism (see Blandford & Ostriker (1978); Bell 
(1978a,b) or Blasi (2013) for a more recent review). In 
particular the excitation of unstable modes by the acceler¬ 
ated particles propagating into the interstellar medium has 
been studied as a possible mechanism to generate strong 
magnetic turbulence predicted by the non-linear diffusive 
shock acceleration theory (see the reviews by Maikov & 
Drury (2001) and Blandford & Eichler (1987)). The study 
of a model problem characterized by two countestreaming 
electron-positron plasma clouds at relativistic energies can 
be relevant to various astrophysical scenarios including 
the fireball model of Gamma Ray Bursts (Piran 2005), 
pulsar wind outflows in Pulsar Wind Nebulae (Blasi & 
Amato 2011), and relativistic jets from Active Galactic 
Nuclei (Begelman et al. 1984). In this context, several au¬ 
thors have studied counterstraming pair plasmas in various 
configurations (see e.g. Hoshino & Shimada (2002); Silva 
et al. (2003); Jaroschek et al. (2005); Ghang et al. (2008); 
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Spitkovsky (2008); Amano & Hoshino (2009); Nishikawa 
et al. (2009); Bret et al. (2013); Liang et al. (2013a,b); 
Lemoine et al. (2014)), including colliding and injected jets 
(as opposed to uniform configurations) which allow the 
generation of collisionless shocks. 

In this paper we examine the instability generated by two 
counter-streaming neutral beams of pair plasmas in the 
ultra-relativistic regime. In particular we address kinetic 
effects, such as particle acceleration, taking place during 
the nonlinear phase of the instability, and we take radiation 
friction (RE) effects into account. It is worth noticing that 
there is a current interest in kinetic simulations of pair 
plasmas with RE included in astrophysics, e.g. for the 
study of anomalous particle acceleration leading to flaring 
in the Grab Nebula (Jaroschek & Hoshino 2009; Gerutti 
et al. 2013). The problem of RE inclusion in the kinetic 
modeling of a relativistic plasma in high electromagnetic 
(EM) fields is also crucial in the context of ultraintense 
laser interaction with matter and plasma (Di Piazza et al. 
2012, and references therein). It is therefore useful to revisit 
classic plasma instabilities in highly relativistic regimes 
possibly dominated by radiation. 

The system composed by two charge-neutral, counter¬ 
streaming pair plasmas is subject to a host of instabilities 
which depend on the orientation of the wavevector with 
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respect to the direction of the beams (for a general review 
see Bret et ah (2010)). The unstable spectrum includes two 
limiting cases: the longitudinal two stream instability (TSI), 
corresponding to an electrostatic mode with flow-aligned 
wavevector, and the transverse filament at ion instability 
(FI), corresponding to an EM mode with wavevector 
perpendicular to the beam direction. The TSI and FI 
are particular cases of the more general instability in 
which the wavevector is oblique to the beam direction 
and the unstable spectrum presents both the electrostatic 
and the EM components. From analytical calculations 
based on first-order perturbation theory (Califano et al. 
(1997); Kazimura et al. (1998); Bret et al. (2004, 2010)), 
the growth rate T in the linear phase for the two-stream 
instability T oc while for the FI T oc where 

7o = (l + (po/'^ec)^) is the initial beam Lorentz factor 
(with po the initial drift momentum). Moreover, these 
calculations show that, when the beams are symmetric, the 
instability is prevalently EM. Thus, in the ultra-relativistic 
regime the transverse FI is expected to dominate the growth 
of the instability, at least before nonlinear effects become 
important. 


We performed EM, fully relativistic particle-in-cell 
(PIC) simulations both in one spatial dimension (ID) 
and in two spatial dimensions (2D) with plane Cartesian 
geometry. In 2D, the simulations can be performed with 
either the counterstreaming beams direction perpendicular 
to the simulation plane (T-mode) or parallel to it (P-mode). 
For the T-mode case, only the growth of transverse modes 
is allowed, while the P-mode allows longitudinal modes 
as well. Thus, in general the dynamics of countestreaming 
instabilities in 2D can be substantially different between T- 
and P-modes (see e.g. Amano & Hoshino (2009) for the case 
of Kelvin-Helmoltz instability in electron-ion plasmas and 
Liang et al. (20I3a,b) for shear instability in pair plasmas) 
so that in principle a three-dimensional (3D) analysis would 
be needed. However, a reliable 3D simulation is often not 
possible because of the huge computational cost, which 
leads to severe limitations in the numerical box size, spatial 
and temporal resolution, and number of particles per cell 
even on a parallel supercomputer. This is particularly 
true for our study where we aim at understanding kinetic 
and particle acceleration effects, which need sufficient 
phase space statistics, i.e. large number of computational 
particles. A similar request holds in order to address the 
effects of the RF force, since the latter is much stronger on 
the highest energy particles in the low-density tail of the 
particle distribution. In addition, the strong coalescence 
of small scale structures in the nonlinear stage eventually 
leads to the formation of structures with size close to the 
numerical box. For these reasons, a “small” 3D simulation 
would excessively suffer from numerical effects at present. 
Therefore in this paper we consider only 2D simulations, 
assuming that a comparison of T- and P-mode simulations 
can give insight into the 3D dynamics. We restrict to a 
configuration of homogeneous, counterstreaming plasmas 
which prevent the formation of shocks, that are not of direct 
interest for this paper. However it should be noticed that 
the nonlinear dynamics and saturation of the instability 


may be different for colliding or injected jets configuration. 

Although as it will be shown below the transverse mode 
is dominant in the early, linear stage leading to the gener¬ 
ation of filaments (which are actually current layers in the 
P-mode), significant differences between the T- and P-mode 
appear in the nonlinear phase. In both cases, the transition 
from the linear to the nonlinear phase is characterized by 
the coalescence of the current filaments, with progressive 
decay of the magnetic field after reaching a peak value at 
the endo of linear phase. Differences between the T- and P- 
mode appear in the nonlinear phase, with the amplitude of 
the magnetic field at peak and at late times being stronger 
for the T-mode. In addition, particle spectra are significantly 
different, with the formation of a spectral peak for the T- 
mode only, while the high energy cut-off is higher for the 
P-mode. Species separation is also different between T- and 
P-modes. The high energy tail of the particle spectrum is 
strongly affected by RF effects, which however do not cause 
substantial modifications in the dynamics of instability and 
in the temporal evolution of fields. 

2 SIMULATION MODEL 
2.1 Numerical set-up 

The initial configuration of our simulations consists of two 
neutral beams of electron-positron pairs which propagate in 
opposite directions (corresponding to pz in the momentum 
space) and fill the entire simulation box. The system is sym¬ 
metric, with the populations of the two beams having the 
same initial density, i.e. 

where tit is the total density, and the same momentum ab¬ 
solute value, i.e. = p^^^ — p ^^2 — Po^ consis¬ 

tently the initial values of charge and current densities and 
of the electric and magnetic fields are zero. A very small 
temperature is introduced to seed the instability. In both 
ID and 2D cases we used periodic boundary conditions. 

We performed simulations with different Lorentz factors 70 
from I to 10^. Here we describe the case with polrUeC = 200 
as it is representative of the most relevant effects observed. 
For different values of po, there are no qualitative changes 
in the dynamics of the instability, the most important dif¬ 
ference being the growth rate of the modes which scales as 
7 q"^^^ (see e.g. Califano et al. (1997); Kazimura et al. (1998); 
Bret et al. (2004)). 

In the ID case, the simulation box is aligned along the x- 
direction and it is divided into 15000 grid cells of equal 
length Ax = O.OI Ap with Xp = c/up the skin depth and 
ujp = [ATve^riT/meY^‘^. Each of the four plasma species 
is represented by Np = 3 x 10® computational particles 
(200 particles per cell). The total simulation time is tsim = 
1000 Tp, where Tp = ^TTjujp^ and the temporal resolution is 
At = O.OITp. 

For the 2D T-mode case, the open-source code PICCANTE 
(Sgattoni et al. 2014; Sgattoni et al. 2015), optimized for 
parallel simulations, has been used. In this case the box had 
2000 X 2000 cells and lengths Lg^x x Lg^y — 100 Ap x 100 Ap 
so Ax = Ap = 0.05 Ap. For each species Np — 2x 10® com¬ 
putational particles and tsim — 200 Tp with At = 0.0325 Tp. 
For the P-mode, simulations performed using the standard 
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Finite Difference Time Domain (FDTD) Maxwell solver al¬ 
gorithm of PICCANTE were strongly affected by numeri¬ 
cal Cherenkov radiation (NCR; for details see Greenwood 
et ah (2004)) due to high-frequency waves which propa¬ 
gate slower than high-energy particles. Thus, for the P-mode 
simulations we set up another PIC code (PICcolino) imple¬ 
menting a spectral Maxwell solver based on the Fast Fourier 
Transform, which is free from NCR. PICcolino was bench- 
marked with PICCANTE in a series of cases where NCR was 
negligible, e.g. in T-mode simulations, showing full agree¬ 
ment in the results. The only noticeable difference was some 
time delay in the early rise of the instability (but with the 
same growth rate) due to a slightly different level of ini¬ 
tial noise. Eor P-mode simulations with PICcolino, the box 
had 1000 X 1000 cells wtih Lg^x x Lg^z = 50Ap x 50 Ap, 
Ax = Az = 0.05 Ap, Ap = 5 X 10^, tsim — 200 Tp and 
At = 0.025 Tp. 

2.2 Radiation friction modeling 

The inclusion of RE in the code is based on the Landau- 
Lifshitz approach (Landau & Lifshitz 1975), with the ap¬ 
proximations and the numerical implementation introduced 
by Tamburini et al. (2010); see also Vranic et al. (2015) for 
a further discussion and comparison to other approaches. 
The radiation friction force which acts on the particles in 
addition to the Lorentz force is 



- (e+^ X b) X B - (^-e) e) , (1) 

where Ve = e^/meC^ ~ 2.8 x 10 “^/im is the classical electron 
radius and A = 27rAp. A third term in the Landau-Lifshitz 
expression has been neglected since it is negligible in all 
situations where use of {rf is appropriate. In the ultra 
relativistic regime the most important contribution in 
Eq. (1) comes from the first term because it is proportional 
to particle Lorentz factor 7 ^ ^ 1. The numerical imple¬ 
mentation in the PIC code is discussed by Tamburini et al. 
(2010). The Compton drag force is neglected. 

In the case without RE inclusion, the equations of the 
PIC code are in an universal dimensionless form with the 
density normalized to nr, time to 1 /cjp, space to c/ojp, and 
fields to mecujp/e. Thus, all the results of a simulation can 
be scaled with respect to a definite value for the density. 
The inclusion of RE breaks such scaling invariance, so 
it is necessary to set a dimensional value for the plasma 
density. We have performed simulations with RE included 
for density values up to tit = 10 ^^cm“^. 

The friction effect of the RE force physically arises from 
the incoherent emission of high-frequency radiation by 
ultra-relativistic electrons and positrons, see Di Piazza 
et al. (2012). Prom a numerical point of view, it is unfea¬ 
sible to perform simulations with a spatial resolution high 
enough to resolve such a small wavelength radiation. Thus, 
it is assumed that such radiation escapes from the system 
without re-interacting with other electrons or positrons, 
and the RE acts as a loss term. Eor density values of 
the order ^ lO^^cm”^, it can be safely assumed that the 
plasma is optically thin to the high-frequency radiation 
(having a typical energy ^ MeV) which mostly contributes 


to radiation losses. In addition, the mean free path for 
Compton scattering of photons is ^ 1 m (as estimated from 
the Klein-Nishina formula), typically much larger than the 
scale length on which the instability sets up (of the order of 
Ip = c/ujp ~ /xm). 


2.3 Symmetry relations 

In a cold four-fluid description, the transverse unstable mode 
exhibits symmetry properties whose violation is a signa¬ 
ture for kinetic and nonlinear effects. Let us indicate the 
density of particles having po > 0 with n% for positrons 
and electrons, respectively, and similarly we define for 
particles having po < 0. We use the same notation for all 
other fluid variables. Eor fields and gradients we use || and 
T to indicate quantities parallel and perpendicular to the 
beams, respectively. Eor the EM transverse unstable mode 
with wavevector k the electric field is parallel to the beams 
(E = (0,0,T^) = E||) while B = {Bx,By,0) = T k 
(we neglect the effect of transverse components of E which 
do not play a role in the linear stage of the instability but 
might be generated due to nonlinear charge separation ef¬ 
fects). The EM field can be thus described via a vector po¬ 
tential A = (0,0, A^) = A|| such that Ey = —^tAy/c and 
B^ = Vx X Ay (we note that V± = {dx ,dy,0) = V). The 
fluid equations can be thus written as 

atni + Vx • (niut«) = 0 , (2) 

AL (7iux.«) = Tie/mec)uf^V± x A,, , (3) 

(liujf;^) = =F(e/mec)L>t7A|| , (4) 

where = {dt + u^V) = {dt + uj +^Vx)- The vector 

potential satisfies the wave equation 

(V^-c-^a?)A|| = -^J|l , (5) 

where 

J|| = + e(ntuj|;^+n+uj|;^) 

-e(n“U|f_^+n“U|f_^) . (6) 

The system is symmetric under the transformation that re¬ 
verses at the same time the charge and the direction of prop¬ 
agation of the populations, which simplifies the description 
of the dynamics. It is also possible to reduce the initial set 
of equations to a system involving only two populations and 
three pairs of dynamical variables, as done in Kazimura et al. 
(1998). The two populations are the sources of the positive 
(Jy") and negative (J^) density current and are identified 
with + and — symbols. The two-fluid variables are defined 
as follows: + nA, ^ = 

uX = u+^ = uX,^, Jj[ = e(n+u+^ - 
and = e(n^Uy"^ — nAu^^). The two-fluid system of 
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Figure 1. ID simulation: the current density Jz as a function 
of position (x) and time (t). The upper plot, panel (a), shows 
the complete evolution over the whole simulation box, while the 
lower plot, panel (b), focuses on the early stage with small-scale 
filaments. The parameters Xp = c/ujp and time Tp = 27r/a;p. 


equations which is obtained from Eqs.(2-4) is 

-I- V± • (n^u*) = 0 , 

( 7 ) 

at(n*u*) -I- Vx • ® uj] = x x Ay) , 

( 8 ) 

2 

[5t-|-Vx •uj]( 7 *jjj=) = + • Vx]A|| . 

(9) 


3 SIMULATION RESULTS 

3.1 One-dimensional simulations 

We first study the FI in ID mostly as a test bed and guid¬ 
ance for multi-dimensional simulations. The ID model has 
the advantages of being directly comparable to analytical 
results, and in particular to check symmetry properties and 
conservation laws. In addition, the ID geometry allows high 
resolution runs and detailed analysis of simulation data. We 
performed several tests changing the number of particles per 
cell and the spatial resolution in order to check the sensitiv¬ 
ity of the results. 

The structure of the current density Jz as a function of 
(x,t) is shown in Fig. 1 (a). The development of the insta¬ 
bility can be divided into three phases: a linear phase for 
t < 50Tp, a transition phase for 50Tp < t < 200Tp and a 
nonlinear, quasi-stationary phase for t > 200Tp. In the lin¬ 
ear phase, modes with a defined wavevector grow exponen¬ 
tially, as we verified by calculating the spatial Fourier trans¬ 



Figure 2. ID simulation: spatial profile of Jz (red line). By (black 
line) and n (total number density, green line) at t = 700Tp. Two 
adjacent maxima or minima of Jz identify a positive or negative 
current filament, respectively. 


form J^[By] = By(kx,t). The numerically obtained growth 
rate for every mode kx agrees well with analytic calculations 
(Kazimura et al. 1998). 

In Fig. 1 (b) a zoom on the structure of Jz during the linear 
and the transition phase is shown. During the exponential 
growth of the perturbations, Jz has a filamentary structure 
with a very small scale length (<C Ap). At t ~ 50 Tp separate 
filaments of opposite current, having a typical scale close to 
the electron skin depth, become distinguishable and start to 
merge. This coalescence characterizes the transition of the 
instability from the linear to the nonlinear quasi-stationary 
regime. For t > 200 Tp the merging phase of the filaments 
finishes and the size of each filament is constant, so that the 
configuration can be described as stationary except for some 
“vibration” which is observable in Fig. 1 (a). 

To understand in more detail the nonlinear regime we con¬ 
sider the spatial profile of Jz, By and n (total number den¬ 
sity), at t — 700Tp, reported in Fig. 2 . A filament with 
positive or negative current is identified by two consecutive 
maxima or minima, respectively. Within each positive fila¬ 
ment, the current density assumes its maximum value near 
the edges. Moving towards the inner region of the filament Jz 
decreases assuming a local minimum at the center, whereas 
the total number density becomes flat-top. The same hap¬ 
pens for negative current filaments. This feature corresponds 
to an anti-correlation between particle density and velocity, 
which will be further discussed below by looking at phase 
space distributions. An oscillatory pattern characterizes also 
the profile of the magnetic field By , which has null points at 
the center of each filament, as it is shown in Fig. 2. 

In the late, quasi-stationary phase the spatial structures of 
Jz, By and n indicate an accumulation of particles within 
the current filaments due to magnetic trapping. In this phase 
the characteristic scale length of the field becomes compa¬ 
rable to the Larmor radius and the density of the magnetic 
energy is of the order of the initial energy density: 

= riT'^QinfieC . ( 10 ) 

OTT 

From Fq. 10 we estimate the Larmor radius as rL,sat = 
y /70 / 2 Ap, which gives the scale length of a filament r/ ~ 
rL,sat ~ 10 Ap, in agreement with the numerical results (see 


JzJenTc) 
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Figure 3. ID simulation: (a) phase space contours {x^pz)\ (b) 
phase space contours (x^px)', (c) effective potential P—(x) = (po + 
a^)^, defined in the eq. (12). Plots (a) and (b) show the behavior 
of the electrons with po/'^eC = 200. All these quantities are 
plotted at t = 800 Tp. The color bar indicates number density. 


Fig. 2). 

Figures 3 (a) and (b) represent, respectively, the projection 
of the phase space on the {x,pz) and the {x,px) planes for 
the electrons with initial positive momentum at t = 800 Tp , 
i.e. described by fluid variables fZ,. The momentum pz is a 
single valued function of the position x^ so that we may also 
speak of pz as a well-defined quantity in fluid equations. In 
the regions of peak density, i.e. in the inner part of each fil¬ 
aments, Pz ~ 0 , which is consistent with the local minimum 
of the current density that peaks at the edge of the filament. 
Outside the filaments there is a small number of electrons 
which have pz ~ 2po = 400. The {x^px) phase space pro¬ 
jection shows a spread along the longitudinal momentum px 
with an approximately Gaussian distribution. 

Consistently with the symmetry properties of the system 
(see Sec. 2 .3), the positrons with initial negative momentum, 
described by /.^, have the same spatial distribution as the 
/G electrons. Thus we may also consider Fig. 3 as being rep¬ 
resentative of the population in the two-fluid description. 
The particles of the f~ population show a pattern analogous 
to Fig. 3 with their spatial distribution in space being com¬ 
plementary to that of the population, i.e. corresponding 
to oppositely directed current filaments. 

In the nonlinear, quasi-st at ionary regime the spatial distri¬ 
bution of particles may be described in terms of an effective 
potential as follows. The conserved canonical momentum is 

= Pz ± ttz, (II) 

where = (e/?TieC^)^z is the dimensionless vector potential 
and zb refers to the sign of the particle charge. At t = 0 we 
have ttz = 0 , so = zbpo for the two beams, respectively. 
The normalized energies (in units of rUeC^) of particles be¬ 


longing to the populations are given by 

= 1 + Px + [TPo + (iz{x)]‘^ = I + + P±{x ). (12) 

The asymptotic state of the system may thus be described 
as a state in which the particles cluster into the minima 
of the effective potential P±{x). Fig. 3 (c) shows P-{x) at 
t — 800 Pp. 

Fig. 4 shows the kinetic energy spectrum of the pop¬ 
ulation for different times, and for both cases in which RF 
is either included or not. As expected from symmetry re¬ 
lations, the spectrum is essentially identical for the other 
three populations. Without RF, the energy spectrum shows 
a sharp a peak at twice the initial kinetic energy, see Fig. 4. 
Correspondingly, we observe a sharp, peaked cut-off at 2po 
in the spectrum of pz (not shown). The peak is strongly 
smoothed in the case with RF, which leads to cooling of 
the plasma by removing particles in the high energy tail, for 
which RF is much stronger due to the ~ 7 ^ scaling. During 
the evolution of the system, RF effects on the particle spec¬ 
tra become more important in the nonlinear phase because 
of the generation of both strong magnetic fields (which lead 
to synchrotron emission) and the acceleration of particles 
to high energy. However, the early development of the in¬ 
stability and the structure and amplitude of the fields at 
saturation are weakly affected by RF. 

The acceleration of particles which double the initial 
value of Pz may be explained as follows. First we notice that 
for each of the two populations (/^ and f~) in the two-fluid 
description of the system, the high-energy particles having 
Pz = =b 2 po are localized outside the filaments where most of 
the particles belonging to the other population are localized, 
as shown in Fig. 3 (a) and Fig. 3 (b). In a given position x 
where the field Ez = Ez{x^ t) acts on a species in such a way 
to reduce its initial momentum po, it necessarily acts on the 
counter-streaming species increasing its initial momentum. 
If a particle belonging to the population falls in a decel¬ 
erating region for the f~ population (i.e. a local minimum 
of the effective potential P-), it gains the same momentum 
Pz that is lost by the particles of the counter-streaming fluid. 

The acceleration mechanism may also be described us¬ 
ing the effective potential, Eq.(I2). For a particle belonging 
to one of the two fluid populations we have 

1 + Px + [4=Po + a{x)]^ = Czp , (13) 

for the f~ and fluids, respectively; are constants. 
Fig. 3 (a-b) shows than for a particle of the f~ population 
Px has a maximum in positions xq where |pF(^o)| = 0, while 
in the same position \pt{xo)\ is maximum and pt{xo) = 0 
for the population. Thus pj assumes its maximum value 
at xo where pF(^o) = 0 i.e. a(xo) = — po- Eq. (13) then 
yields C+ = I + 4po. Due to the symmetry of the system, 
the vector potential a{x) assumes the same values for its 
maxima and minima as a function of x, so there will be 
another point Xq where a(xo) = po, which using Eq. (II) 
yields the maximum value of the momentum |pj (xo)| = 2po. 
Eor symmetry reasons we also obtain C- = I + 4po and a 
maximum of 2po for |pF | • Thus the maximum energy of the 
particles is 

£’„ax = (l+4pg)'/". (14) 

The particles with maximum pz are in positions where 
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Figure 4. ID simulation: kinetic energy spectrum of the pop¬ 
ulation at different times, for simulations with and without RF. 
In the case without RF (shown at t = 400 Tp) a spectral peak 
appears at twice the initial beam energy ~ 27 omeC^, and the 
spectrum is almost unchanged at late times. For the case with 
RF, the spectral peak is smoothed out at t = 400 Tp and at later 
times {t = 1000 Tp) the high energy tail is “washed” out because 
of radiative losses. 

has a maximum or minimum. Thus, dpzjdx — -^dazjdx = 0 
where \pz\ = 2po, i-e. the particles all gain the same mo¬ 
mentum to first order in their distance from the maximum 
of ja^l, which explains the peak at the cut-off in the energy 
spectrum (Fig. 4) as the formation of a spectral caustic. 

3.2 Two-dimensional simulations 

In this section we present the 2D simulations, comparing 
the results of T-mode and P-mode geometry. As already 
mentioned in Sec.2.1, two different codes have been used 
for numerical reasons. The rise of the instability is shifted 
in time between T- and P-mode simulations because of the 
different noise level in the two codes, although the growth 
rate is identical in benchmark cases. Thus, to make com¬ 
parisons at the same physical time, the simulation time has 
been shifted in order that the instant at which the mag¬ 
netic energy reaches its peak (marking the end of the linear 
growth stage) coincides for the two cases. 

Fig.5 shows the distribution of Jz (the current density par¬ 
allel to the beams direction) for the P- and T—plane cases, 
at three different times. In order to also represent the mag¬ 
netic field distribution. Fig.6 shows Bz in the P-plane and 
{Bx + in the T-plane at the same times of Fig.5. 

In the early stage of the simulations (t < 50), the P-case 
shows parallel current filaments which are elongated in the 
beam direction z and have almost the same width as the 
transverse structures in the T-case. This confirms that in 
the linear stage the most unstable wavevector is along the 
direction of the beams, i.e. the FI is of transverse nature. 
During this early stage, the amplitude of the field grows 
exponentially. At saturation (t ~ 55), the beam energy con¬ 
verted into magnetic energy for the T-case is nearly two 
times the value for the P-case (Fig.7), in agreement with an 
approximate energy equipartition. 

At later times (t > 55) merging of small-scale filaments is 
observed in both the T- and P-cases, eventually leading at 
long times (t > 100) to the formation of structures with a 
size close to that of the numerical box for both cases. How¬ 
ever, significant differences are apparent between the T- and 
P-cases. 


For the T-case, the current distribution across a large-scale 
2D structure is similar to that observed in ID: the current 
peaks near the boundary of the island (which corresponds 
to the “horned” ID profiles in Fig.2) and has much weaker 
values well inside the island; locally, small scale filaments 
where the current changes sign are also observed, which are 
caused by the different orientation of the wave-vector with 
respect to the initial direction of the beams. The magnetic 
field is strongly localized along the boundary of the current 
structure, i.e. along null lines of Jz- The spatial correlation 
between the density of each species and the fields is also 
similar to that observed in ID. The current and density dis¬ 
tributions during the non-linear phase are similar to those 
observed as asymptotic numerical solutions of 2D Navier- 
Stokes and magnetohydrodynamic equations (Hossain et al. 
1983). The large scale structures of the magnetic field evolve 
slowly both in the shape and in the amplitude of the field 
for t > 100. 

The distributions of total kinetic energy and pz for the T- 
case are shown in Fig. 8 (frame (5) and {d)). As in the ID 
case, a peak at the upper cut-off pz — 2pQ — 400meC forms 
(see Fig. 8 (d)), while the spectral peak in the energy dis¬ 
tribution (Fig. 8 (b)) disappears, because the energy in the 
tail of the distribution is “smoothed” out over the additional 
degree of freedom. Hence we can expect that the spectrum 
would be smeared out in a 3D case. While the inclusion of 
RF is found not to change the growth and development of 
the filaments significantly, it has a major impact on the high 
energy tail of the spectrum, reducing the cut-off by a factor 
of ~ 2, similarly to what observed in ID. The amount of 
energy lost to radiation exceeds 30% at the end of the sim¬ 
ulation (Fig.9). 

In the P-case, bending and tearing of filaments during the 
merging stage is observed. This leads to the generation of 
structures which are not homogeneous along the direction 
of the beams, i.e. to a spectrum of modes with kz % 0. The 
latter processes can not be simply viewed as the growth 
of an unstable longitudinal mode: a Fourier analysis high¬ 
lights a broad spectrum in kz at late times. Correspondingly, 
electrostatic fields are generated leading to breaking of the 
symmetry properties of the system for purely transverse EM 
perturbations. 

In the P-case the large scale structures of the magnetic field 
are less regular than in the T-case, showing a small-scale 
irregular structure at t = 200. The decay of the magnetic 
field is much more pronounced with the magnetic energy 
becoming of the order of the electrostatic energy at the end 
of the simulation (t — 200). This behavior is likely to be 
due to electrostatic fields causing heating of electrons and 
positrons in the simulations plane. The energy spectrum in 
the P-case becomes broader than in the T-case, with an 
higher energy cut-off. No narrow peak is observed at 2po in 
the P-case, confirming that peak formation is related to the 
conservation of canonical momentum along the direction of 
the beams in the T-case. 

Finally, we discuss the importance of radiative losses due to 
the inclusion of RF for different geometry and density. Fig¬ 
ure 9 reports the time evolution of the total energy (particle 
energy plus field energy) with respect to the initial kinetic 
energy of the beams, for different simulations with RF in¬ 
cluded. Radiative losses are higher for the T-mode than for 
the P-mode, consistently with the higher fields generated in 
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Figure 5. Spatial distribution of Jz in 2D simulations for the P-mode (lower row) and T-mode (upper row) cases, at four different times 
(simulation times have been shifted in order than the instant of peak magnetic energy coincides for the two cases, see text for details). 
Only a quarter of the whole box is plotted for the T-mode simulation. 


the T-mode case. The effect of RF is also stronger for higher 
densities nr , which is also consistent with the magnetic field 
at saturation being proportional to At densities of the 
order of 10 ^^cm“^, radiative losses reach a few per cent of 
the initial energy at the end of the simulation. 

Although for very high density there is a major loss of en¬ 
ergy due to the RF effect, the instability dynamics is not 
strongly changed with respect to the case without RF. The 
system organizes itself in filamentary structures for the cur¬ 
rent density which have almost the same size and features 
of the filaments obtained in the non-RF simulation. This 
behavior can be simply understood by noticing that the 
EM fields have to grow in order for RF to be important, 
so that the RF plays little role before the saturation phase. 
Moreover, in the ultra-relativistic case the dominant term of 
RF force (see Landau & Lifshitz (1975)) is oc 7 ^. Thus the 
RF contribution is strongly increased by the acceleration of 
some particles to higher energy, which is maximized at the 
instability saturation stage. This is consistent with RF ef¬ 
fects being more evident in the particle spectra, as shown in 
Fig. 8 . 

4 CONCLUSIONS 

In this work we have studied the evolution of the filamen¬ 
tation instability produced by two counter-streaming pair 


plasmas using PIC simulations in ID and in 2D for both T- 
and P-modes, with and without radiation friction effects. 
The saturation level of the instability and the particle spec¬ 
tra are significantly different between T- and P-modes. In 
the T-mode case, the magnetic field at saturation is stronger 
and has a slower decay in time after reaching its maximum 
value; the particle spectrum shows the formation of a spec¬ 
tral peak at cut-off for which a simple theory has been pre¬ 
sented. In the P-mode case, the magnetic field has a lower 
maximum value and has a faster decay, so that the mag¬ 
netic energy becomes comparable to the electrostatic energy 
at the end of the simulations; the energy spectra show no 
peak but a higher energy cut-off. Radiation friction effects 
have been found to be strong only for relatively high density 
(~ I 0 ^°cm“^) and to modify strongly the particle spectra, 
cooling down the distribution functions and removing the 
highest energy particles, while the instability development 
is weakly affected. 
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